library(lavaan)
library(foreign)
setwd("~/Replication Code/")

##---Loading in data ----
anes9297 <- read.dta("./data/anes9297_subset.dta")

ccap08 <- read.dta("./data/ccap08_subset.dta")

ccap12 <- read.dta("./data/ccap12_subset.dta")

vsg.dat <- read.dta("./data/vsg_subset.dta")

ccap16 <- read.dta("./data/ccap16_subset.dta")

FIT <- c("chisq","df","cfi", "tli", "srmr", "rmsea", "rmsea.ci.lower", "rmsea.ci.upper")

#----------------------------------------------------------------------------------------------------#
# ANES 1992-1994
#----------------------------------------------------------------------------------------------------#
m.anes.sem <- '
# measurement model
rr92 =~ NA*rr_irish_92 + v1*rr_irish_92 + v2*rr_blklsds_92 + v3*rr_tryhrd_92 + v4*rr_slvry_92
rr94 =~ NA*rr_irish_94 + v1*rr_irish_94 + v2*rr_blklsds_94 + v3*rr_tryhrd_94 + v4*rr_slvry_94

pid92 =~ pid92_7scaled
pid94 =~ pid94_7scaled
# regressions
rr94 ~ a*pid92 + rr92 
pid94 ~ pid92 + b*rr92
# coef differences?
coef.dif := a-b

# setting variances
rr92 ~~ 1*rr92
rr94 ~~ 1*rr94

# wording item error variances
#rr_irish_92 ~~ rr_tryhrd_92
rr_irish_94 ~~ rr_tryhrd_94

rr_blklsds_92 ~~ rr_slvry_92
rr_blklsds_94 ~~ rr_slvry_94

# temporal item error variances
#rr_irish_92 ~~ rr_irish_94
rr_tryhrd_92 ~~ rr_tryhrd_94
rr_blklsds_92 ~~ rr_blklsds_94
rr_slvry_92 ~~ rr_slvry_94

# constraining covariances
rr92 ~~ cov*pid92
rr94 ~~ cov*pid94
'

dat1 <- subset(anes9297, white_int92 == white_int94,
               select = c("pid92_7scaled", "pid94_7scaled",
                          "V940005",
                          "rr_irish_92", "rr_blklsds_92", "rr_tryhrd_92", "rr_slvry_92",
                          "rr_irish_94", "rr_blklsds_94", "rr_tryhrd_94", "rr_slvry_94"))
dat1 <- na.omit(dat1)
fit.anes9294.sem <- sem(m.anes.sem, data = dat1, mimic = "Mplus")

fit.anes9294 <- fit.anes9294.sem
summary(fit.anes9294, standardized = T)

fit_main <- fitMeasures(fit.anes9294)[FIT]
fit_main

#----------------------------------------------------------------------------------------------------#
# 2008 CCAP
#----------------------------------------------------------------------------------------------------#
m.ccap08.sem <- '
# measurement model (incorporate first item twice because parse works on single options)
rrM =~ NA*M_specfav + v1*M_specfav + v2*M_deserve + v3*M_tryhard + v4*M_pdisc
rrO =~ NA*O_specfav + v1*O_specfav + v2*O_deserve + v3*O_tryhard + v4*O_pdisc
pidM =~ pid7r_m_sc
pidO =~ pid7r_o_sc
# regressions
rrO ~ a*pidM + rrM
pidO ~ pidM + b*rrM
# coef differences?
coef.dif := a-b

# setting variances
rrM ~~ 1*rrM
rrO ~~ 1*rrO

# wording item error variances
M_specfav ~~ M_tryhard
O_specfav ~~ O_tryhard
M_deserve ~~ M_pdisc
O_deserve ~~ O_pdisc

# temporal item error variances
M_specfav ~~ O_specfav
M_tryhard ~~ O_tryhard
M_deserve ~~ O_deserve
M_pdisc ~~ O_pdisc

# constraining covariances
rrM ~~ cov*pidM
rrO ~~ cov*pidO
'

dat1 <- subset(ccap08,
               select = c("pid7r_m_sc", "pid7r_o_sc",
                          "weight",
                          "M_specfav", "M_deserve", "M_tryhard", "M_pdisc",
                          "O_specfav", "O_deserve", "O_tryhard", "O_pdisc"))
dat1 <- na.omit(dat1)
fit.ccap08.sem <- sem(m.ccap08.sem, data = dat1, mimic = "Mplus")
fit.ccap08 <- fit.ccap08.sem
summary(fit.ccap08, standardized = T)

fit_main <- fitMeasures(fit.ccap08)[FIT]
fit_main

#----------------------------------------------------------------------------------------------------#
# CCAP 2012
#----------------------------------------------------------------------------------------------------#
## Period 1
m.ccap12p1.sem <- '
# measurement model (incorporate first item twice because parse works on single options)
rr_1 =~ NA*pp_specfav + v1*pp_specfav + v2*pp_deserve + v3*pp_tryhard + v4*pp_pdisc
rr_2 =~ NA*specfav + v1*specfav + v2*deserve + v3*tryhard + v4*pdisc
pid_1 =~ pp_pid7r_sc
pid_2 =~ pid7_sc
# regressions
rr_2 ~ a*pid_1 + rr_1 
pid_2 ~ pid_1 + b*rr_1
# coef differences?
coef.dif := a-b

# setting variances
rr_1 ~~ 1*rr_1
rr_2 ~~ 1*rr_2

# wording item error variances
pp_specfav ~~ pp_tryhard
specfav ~~ tryhard
pp_deserve ~~ pp_pdisc
deserve ~~ pdisc

# temporal item error variances
pp_specfav ~~ specfav
pp_tryhard ~~ tryhard
pp_deserve ~~ deserve
pp_pdisc ~~ pdisc

# constraining covariances
rr_1 ~~ cov*pid_1
rr_2 ~~ cov*pid_2
'

dat1 <- subset(ccap12, p2 == 1, 
               select = c("pp_pid7r_sc", "pid7_sc",
                          "weight",
                          "pp_specfav", "pp_deserve", "pp_tryhard", "pp_pdisc",
                          "specfav", "deserve", "tryhard", "pdisc"))
dat1 <- na.omit(dat1)
fit.ccap12p1.sem <- sem(m.ccap12p1.sem, data = dat1, mimic = "Mplus")
fit.ccap12p1 <- fit.ccap12p1.sem
summary(fit.ccap12p1, standardized = T)

fit_main <- fitMeasures(fit.ccap12p1)[FIT]
fit_main



####### Period 2
m.ccap12p2.sem <- '
# measurement model (incorporate first item twice because parse works on single options)
rr_1 =~ NA*pp_specfav + v1*pp_specfav + v2*pp_deserve + v3*pp_tryhard + v4*pp_pdisc
rr_2 =~ NA*specfav + v1*specfav + v2*deserve + v3*tryhard + v4*pdisc
pid_1 =~ pp_pid7r_sc
pid_2 =~ pid7_sc
# regressions
rr_2 ~ a*pid_1 + rr_1 
pid_2 ~ pid_1 + b*rr_1
# coef differences?
coef.dif := a-b

# setting variances
rr_1 ~~ 1*rr_1
rr_2 ~~ 1*rr_2

# wording item error variances
pp_specfav ~~ pp_tryhard
specfav ~~ tryhard
pp_deserve ~~ pp_pdisc
deserve ~~ pdisc

# temporal item error variances
pp_specfav ~~ specfav
pp_tryhard ~~ tryhard
pp_deserve ~~ deserve
pp_pdisc ~~ pdisc

# constraining covariances
rr_1 ~~ cov*pid_1
rr_2 ~~ cov*pid_2
'

dat1 <- subset(ccap12, p3 == 1, 
               select = c("pp_pid7r_sc", "pid7_sc",
                          "weight",
                          "pp_specfav", "pp_deserve", "pp_tryhard", "pp_pdisc",
                          "specfav", "deserve", "tryhard", "pdisc"))
dat1 <- na.omit(dat1)
fit.ccap12p2.sem <- sem(m.ccap12p2.sem, data = dat1, mimic = "Mplus")

fit.ccap12p2 <- fit.ccap12p2.sem
summary(fit.ccap12p2, standardized = T)

fit_main <- fitMeasures(fit.ccap12p2)[FIT]
fit_main

#------------------------------------------------------------------------------#
# 2012-2016 Voter Survey
#----------------------------------------------------------------------------------------------------#
m.vsg.sem <- '
# measurement model
rr12 =~ NA*rr_specfav_12 + v1*rr_specfav_12 + v2*rr_deserve_12 + v3*rr_thard_12 + v4*rr_pdisc_12
rr16 =~ NA*rr_specfav_16 + v1*rr_specfav_16 + v2*rr_deserve_16 + v3*rr_thard_16 + v4*rr_pdisc_16
pid12 =~ pid7_sc_12
pid16 =~ pid7_sc_16

# regressions
rr16 ~ a*pid12 + rr12
pid16 ~ pid12 + b*rr12
# coef differences?
coef.dif := a-b

# setting variances
rr12 ~~ 1*rr12
rr16 ~~ 1*rr16

# wording item error variances
rr_specfav_12 ~~ rr_thard_12
rr_specfav_16 ~~ rr_thard_16
rr_deserve_12 ~~ rr_pdisc_12
rr_deserve_16 ~~ rr_pdisc_16

# temporal item error variances
rr_specfav_12 ~~ rr_specfav_16
rr_thard_12 ~~ rr_thard_16
rr_deserve_12 ~~ rr_deserve_16
rr_pdisc_12 ~~ rr_pdisc_16

# constraining covariances
rr12 ~~ cov*pid12
rr16 ~~ cov*pid16
'


dat1 <- subset(vsg.dat,
               select = c("pid7_sc_12", "pid7_sc_16",
                          "weight",
                          "rr_deserve_12", "rr_specfav_12", "rr_thard_12", "rr_pdisc_12",
                          "rr_deserve_16", "rr_specfav_16", "rr_thard_16", "rr_pdisc_16"))
dat1 <- na.omit(dat1)

fit.vsg.sem <- sem(m.vsg.sem, data = dat1, mimic = "Mplus")

fit.vsg <- fit.vsg.sem
summary(fit.vsg, standardized = T)

fit_main <- fitMeasures(fit.vsg)[FIT]
fit_main


#----------------------------------------------------------------------------------------------------#
# 2016 CCAP
#----------------------------------------------------------------------------------------------------#
m.ccap16.sem <- '
# measurement model
rrB =~ NA*b_ldes + v2*b_ldes + v3*b_thard + v4*b_pdisc
rrP =~ NA*p_ldes + v2*p_ldes + v3*p_thard + v4*p_pdisc
pidB =~ b_pid7_sc
pidP =~ p_pid7_sc
# regressions
rrP ~ a*pidB + rrB
pidP ~ pidB + b*rrB
# coef differences?
coef.dif := a-b

# setting variances
rrB ~~ 1*rrB
rrP ~~ 1*rrP

# wording item error variances
b_ldes ~~ b_pdisc
p_ldes ~~ p_pdisc

# temporal item error variances
b_thard ~~ p_thard
b_ldes ~~ p_ldes
b_pdisc ~~ p_pdisc

# constraining covariances
rrB ~~ cov*pidB
rrP ~~ cov*pidP
'

dat1 <- subset(ccap16,
               select = c("b_pid7_sc", "p_pid7_sc",
                          "weight_post",
                          "b_ldes", "b_thard", "b_pdisc",
                          "p_ldes", "p_thard", "p_pdisc"))
dat1 <- na.omit(dat1)
fit.ccap16.sem <- sem(m.ccap16.sem, data = dat1, mimic = "Mplus")

fit.ccap16 <- fit.ccap16.sem
summary(fit.ccap16, standardized = T)
fit_main <- fitMeasures(fit.ccap16)[FIT]
fit_main


#----------------------------------------------------------------------------------------------------#
# Tables
#----------------------------------------------------------------------------------------------------#
## Measurement Model Results
sem.models <- c("fit.anes9294.sem", "fit.ccap08.sem", "fit.ccap12p1.sem", "fit.ccap12p2.sem",
                "fit.vsg.sem", "fit.ccap16.sem")

items <- 17
OUT.MEASURE <- array(NA, c(items, length(sem.models)))
rownames(OUT.MEASURE) <- c("Special Favors", "se sf",
                           "Deserve Less", "se dl",
                           "Try Hard", "se th",
                           "Past Discrimination", "se pd",
                           "Partisanship", "se pid",
                           "Chi-Square", "DF", "CFI", "TLI", "SRMR", "RMSEA [90% CI]", "N")

OUT.MEASURE <- as.data.frame(OUT.MEASURE)


for(i in 1:length(sem.models)){
  df <- get(sem.models[i])
  est <- parameterEstimates(df)
  # Factor Loadings
  if("v1" %in% est$label){
    OUT.MEASURE["Special Favors", i] <- round(est$est[1], 3)
    OUT.MEASURE["se sf", i] <- paste0("(", round(est$se[1], 3), ")")
    OUT.MEASURE["Deserve Less", i] <- round(est$est[2], 3)
    OUT.MEASURE["se dl", i] <- paste0("(", round(est$se[2], 3), ")")
    OUT.MEASURE["Try Hard", i] <- round(est$est[3], 3)
    OUT.MEASURE["se th", i] <- paste0("(", round(est$se[3], 3), ")")
    OUT.MEASURE["Past Discrimination", i] <- round(est$est[4], 3)
    OUT.MEASURE["se pd", i] <- paste0("(", round(est$se[4], 3), ")")
  }
  else{
    OUT.MEASURE["Deserve Less", i] <- round(est$est[2], 3)
    OUT.MEASURE["se dl", i] <- paste0("(", round(est$se[2], 3), ")")
    OUT.MEASURE["Try Hard", i] <- round(est$est[3], 3)
    OUT.MEASURE["se th", i] <- paste0("(", round(est$se[3], 3), ")")
    OUT.MEASURE["Past Discrimination", i] <- round(est$est[4], 3)
    OUT.MEASURE["se pd", i] <- paste0("(", round(est$se[4], 3), ")")
  }
  OUT.MEASURE["Partisanship", i] <- 1
  OUT.MEASURE["se pid", i] <- "(---)"
  
  ## Folding in fit measures
  OUT.MEASURE["Chi-Square", i] <- round(fitMeasures(df)["chisq"], 3)
  OUT.MEASURE["DF", i] <- round(fitMeasures(df)["df"], 3)
  OUT.MEASURE["CFI", i] <- round(fitMeasures(df)["cfi"], 3)
  OUT.MEASURE["TLI", i] <- round(fitMeasures(df)["tli"], 3)
  OUT.MEASURE["SRMR", i] <- round(fitMeasures(df)["srmr"], 3)
  OUT.MEASURE["RMSEA [90% CI]", i] <- paste0(round(fitMeasures(df)["rmsea"], 3),
                                             " [",
                                             round(fitMeasures(df)["rmsea.ci.lower"], 3),
                                             ", ",
                                             round(fitMeasures(df)["rmsea.ci.upper"], 3),
                                             "]")
  OUT.MEASURE["N", i] <- nobs(df)
  
}
tab <- xtable::xtable(OUT.MEASURE, 
                      label = c("measure_mods"), 
                      align = c("l", rep("c", ncol(OUT.MEASURE))),
                      caption = c("Measurement Model Results"))
print(tab)

#### STRUCTURAL RESULTS

# Creating Placehold Data Frames for Tables
OUT.stab <- array(NA, c(4, length(sem.models)*2))
OUT.clag <- array(NA, c(4, length(sem.models)*2))
rownames(OUT.stab) <- c("Partisanship", "pid se", "Racial Resentment", "rr se")
OUT.stab <- as.data.frame(OUT.stab)
rownames(OUT.clag) <- c("Partisanship", "pid se", "Racial Resentment", "rr se")
colnames(OUT.stab) <- c("ANES 1992-1994", "", "CCAP 2008", "",
                        "CCAP 2012: March", "", "CCAP 2012: August", "",
                        "2012-2016 VOTER Survey", "", "CCAP 2016", "")
OUT.stab <- as.data.frame(OUT.stab)
rownames(OUT.clag) <- c("Partisanship", "pid se", "Racial Resentment", "rr se")
colnames(OUT.clag) <- c("ANES 1992-1994", "", "CCAP 2008", "",
                        "CCAP 2012: March", "", "CCAP 2012: August", "",
                        "2012-2016 VOTER Survey", "", "CCAP 2016", "")
OUT.clag <- as.data.frame(OUT.clag)

# Exporting Results
k <- 1 #indexer
for(i in 1:length(sem.models)){
  df <- get(sem.models[i])
  est <- parameterEstimates(df)
  RR.vars <- unique(grep("rr", est$lhs, value = T))[which(nchar(unique(grep("rr", est$lhs, value = T))) < 8)]
  pid.vars <- unique(grep("pid", est$lhs, value = T))[which(nchar(unique(grep("pid", est$lhs, value = T))) < 8)]
  ## Stablities
  OUT.stab["Partisanship", k] <- round(est$est[which(est$lhs == pid.vars[2] & est$rhs == pid.vars[1] & est$op == "~")], 3)
  OUT.stab["pid se", k] <- paste0("(", round(est$se[which(est$lhs == pid.vars[2] & est$rhs == pid.vars[1] & est$op == "~")], 3), ")")
  OUT.stab["Racial Resentment", k] <- round(est$est[which(est$lhs == RR.vars[2] & est$rhs == RR.vars[1]  & est$op == "~")], 3)
  OUT.stab["rr se", k] <- paste0("(", round(est$se[which(est$lhs == RR.vars[2] & est$rhs == RR.vars[1]  & est$op == "~")], 3), ")")
  ## Cross-lagged effects
  OUT.clag["Partisanship", k] <- round(est$est[which(est$lhs == RR.vars[2] & est$rhs == pid.vars[1] & est$op == "~")], 3)
  OUT.clag["pid se", k] <- paste0("(", round(est$se[which(est$lhs == RR.vars[2] & est$rhs == pid.vars[1] & est$op == "~")], 3), ")")
  OUT.clag["Racial Resentment", k] <- round(est$est[which(est$lhs == pid.vars[2] & est$rhs == RR.vars[1]  & est$op == "~")], 3)
  OUT.clag["rr se", k] <- paste0("(", round(est$se[which(est$lhs == pid.vars[2] & est$rhs == RR.vars[1]  & est$op == "~")], 3), ")") 
  
  ### Completely Standardized Solution
  est.std <- standardizedSolution(df)
  ## Stablities
  OUT.stab["Partisanship", (k+1)] <- paste0("\\textit{", round(est.std$est.std[which(est.std$lhs == pid.vars[2] & est.std$rhs == pid.vars[1] & est.std$op == "~")], 3),"}")
  OUT.stab["pid se", (k+1)] <- paste0("\\textit{(", round(est.std$se[which(est.std$lhs == pid.vars[2] & est.std$rhs == pid.vars[1] & est.std$op == "~")], 3), ")}")
  OUT.stab["Racial Resentment", (k+1)] <- paste0("\\textit{", round(est.std$est.std[which(est.std$lhs == RR.vars[2] & est.std$rhs == RR.vars[1]  & est.std$op == "~")], 3),"}")
  OUT.stab["rr se", (k+1)] <- paste0("\\textit{(", round(est.std$se[which(est.std$lhs == RR.vars[2] & est.std$rhs == RR.vars[1]  & est.std$op == "~")], 3), ")}")
  ## Cross-lagged effects
  OUT.clag["Partisanship", (k+1)] <- paste0("\\textit{", round(est.std$est.std[which(est.std$lhs == RR.vars[2] & est.std$rhs == pid.vars[1] & est.std$op == "~")], 3),"}")
  OUT.clag["pid se", (k+1)] <- paste0("\\textit{(", round(est.std$se[which(est.std$lhs == RR.vars[2] & est.std$rhs == pid.vars[1] & est.std$op == "~")], 3), ")}")
  OUT.clag["Racial Resentment", (k+1)] <- paste0("\\textit{", round(est.std$est.std[which(est.std$lhs == pid.vars[2] & est.std$rhs == RR.vars[1]  & est.std$op == "~")], 3),"}")
  OUT.clag["rr se", (k+1)] <- paste0("\\textit{(", round(est.std$se[which(est.std$lhs == pid.vars[2] & est.std$rhs == RR.vars[1]  & est.std$op == "~")], 3), ")}")
  
  # upping the indexer
  k <- k + 2
}


stab.tab <- xtable::xtable(OUT.stab, 
                           label = c("stabilities"), 
                           align = c("l", rep("c", ncol(OUT.stab))),
                           caption = c("Stability Coefficients for Partisanship and Racial Resentment"))
print(stab.tab)

clag.tab <- xtable::xtable(OUT.clag, 
                           label = c("causal_effects"), 
                           align = c("l", rep("c", ncol(OUT.clag))),
                           caption = c("Cross-Lagged Effects of Partisanship and Racial Resentment"))
print(clag.tab)
